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ABSTRACT 

To be presented is a study of the secular evolution of a spherical stellar system with 
a central star-accreting black hole (BH) using the anisotropic gaseous model. This 
method solves numerically moment equations of the full Fokkcr-Planck equation, with 
Boltzmann-Vlasov terms on the left- and collisional terms on the right-hand sides. We 
study the growth of the central BH due to star accretion at its tidal radius and the 
feedback of this process on to the core collapse as well as the post-collapse evolution 
of the surrounding stellar cluster in a self-consistent manner. Diffusion in velocity 
space into the loss-cone is approximated by a simple model. The results show that the 
self-regulated growth of the BH reaches a certain fraction of the total mass cluster 
and agree with other methods. Our approach is much faster than competing ones 
(Monte Carlo, iV-body) and provides detailed informations about the time and space 
dependent evolution of all relevant properties of the system. In this work we present 
the method and study simple models (equal stellar masses, no stellar evolution or 
collisions). Nonetheless, a generalisation to include such effects is conceptually simple 
and under way. 

Key words: anisotropy - star clusters - stellar dynamics - black holes 



1 INTRODUCTION 

The quest for the source of the luminosities of L w 10 12 Lq 
produced on very small scales, jets and other properties of 
quasars and other types of active galactic nuclei (AGN) led 
in the 60's and 70's to a thorough research that hint to 
the inkling of "super-massive central objects" harboured at 
their centres. These were suggested to be the main source of 
such characteristics <Lvnden-Belil967tlLvnden-Bell fc Reesl 
ll97ll:lHillJl975h . Two years later. lLvnden-Belll96Sl showed 
that the release of gravitational binding energy by stellar 
accretion on to a SMBH could be the primary powerhouse 
of an AGN iLvnden-Bellll969lh Nowadays, t his idea see ms 
to be ensconced in extra galactic astronomy dReesl |l984'). 

A direct consequence of the paradigm of SMBHs at the 
centre of ancient galaxies to explain the energy emitted by 
quasars is that relic SMBHs s hould inha bit at least a frac- 
tion of present-day galaxies jReeslll99 (J). This conclusion 
was first made quantitative b^ Tsoltanl <I1982|) and has re- 
cently be revisi ted in more deta i l and in the light of recent 
observations bv lYu fc Tremainel i2002t) . 

In the last decade, observational evidences have been 



accumulating that strongly suggest that massive BHs are 
indeed present at the centre of most galaxies, with a signifi- 
cant spheroidal component. Mostly thanks to the HST, the 
kinematics in present-day universe of gas or stars has been 
measured in the central parts of tens of nearby galaxies. In 
almost all cases , proper modelling of the measured mo- 
tions requires the presence of a centra l compact dark object 
with a mass of a few 10 6 to 10 9 Mr^ ifFerrarese etld] | 200l| : 
iGebhardt et al.ll2002l: iPinknev et"al]|2003l : lKormendvll2005 , 



and references therein). Note, however, that the conclusion 
that such an object is indeed a BH rather than a cluster 
of smaller dark objects (like neutron stars, brown dwarves 
etc) has only been reached for two galaxies. The first one 
is the Milky Way itself at the centre of which the case for 
a 3-4 x 10 6 M Q MBH has been clinched, mostly through 
ground-based I R observations o f the fast orbital m otions 
of a few stars jGhez et alJl2003t ISchodel et alj|2003lh The 
second case is NGC4258, which passes a central Keplerian 
gaseous disc with H2O MASER strong sources allowing high 
resolution VLBI observations down to 0.16 pc of the cen- 
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tre llMivoshi et alll99llHerrnstein et alll999l;lMoran et alJ 
Il999h . 

In any case, it is nowadays largely accepted that the cen- 
tral dark object required to explain kinematics data in local 
active and non-active galaxies is an MBH. The large number 
of galaxies surveyed has allowed to study the demographics 
of the MBHs and, in particular, look for correlations with 
properties of the host galaxy. The most remarkable ones are 
the fact that the MBH has a mass which is roughly about 
0.1% of the stellar mass of the spheroidal component of the 
galaxy and that the mass of the BH, Mbh, correlates even 
more tightly with the velocity of this component. These facts 
certainly strike a close link between the formation of the 
galaxy and the massive object harboured at its centre. If one 
extends these relations to smaller stellar systems, one could 
expect that globular clusters host so-called intermediate- 
mass black holes, i.e. BHs whose mass is in the range of 
10 2 -10 4 M Q . After having been suggested in the 70's to ex- 
plain the x-ray sources observed in globulars clusters, later 
discovered to be stellar-mass binaries, this possibility has re- 
cently be revived by two lines of observations. First IMBHs 
may explain the ultra-luminous x-ray sources (ULXs) that 
are present in regions of strong stellar formation in inter- 
acting galaxies and hence suggesting a link with young "su- 
per stellar clusters" (SSC), although ULXs are typically not 
found at the centre of luminous SSCs. On IMBHs and their 
possible link to ULXs, see the review bv iMiller fc Colbertl 
i2003l) . Second, recent HST observations of the stellar kine- 
matics at the centre of M15 around the Milky Way and Gl 
around M31 have been interprete d as indications of the pres- 
ence of an IMBH in both clusters (Ivan der Marel et al .120021: 
iGerssen etalll20ul 120031: iGebhardt et al.ll2002f) . However, 
in the case of M15, the mass of the point masses required by 
the observations is compatible with zero and TV-body mod- 
els have been made of both clusters that la ck a central IMBH 
but ar e compatible with the observations feaumgardt et all 
l2003al lbTl . We note that scenarios have been proposed that 
would quite naturally explain the formation of an IMBH at 
the centre of a stellar cluster, through run-away stellar colli- 
sions, provided that the relaxation time is short enough and 



that very massive stars 


10 2 Mo < M* < 10 4 M(7,) evolve into 


IMBH fcbisuzaki et al. 


1200 ll IPorteeies Zwart & McMillan! 


20021: iRasio et alj|2003l: 


Giirkan et alll2004). 



The theoretical study of the structure and evolution of 
a stellar cluster (galactic nucleus or globular cluster) har- 
bouring a central MBH started 30 years ago. However, due 
to the complex nature of the problem which includes many 
physical processes and span a huge range of time and length 
scales, our understanding of such systems is still incomplete 
and, probably, subjected to revision. As in many fields of 
astrophysics, analytical computations can only been applied 
to highly idealised situations and only a very limited variety 
of numerical methods have been developed so far that can 
tackle this problem. 

In the present paper, we introduce a simulation method 
to follow the joint evolution of a spherical star cluster with 
a central BH making feasible anisotropy. This procedure is 
based on moments of the Boltzmann equation with relax- 
ation. The cluster is modelled like a self-gravitating, con- 
d ucting gas sphere, accor ding to the methods presente d 
in lLouis fc Spurzeml dl99lT) and iGiersz fc Spurzeml lll994D. 
These models improved earlier gas models of iBettwieserl 



lll983ft and lHeggiel lll984h . Much as the structure of the nu- 
merical method is for the sake of computational efficiency 
on account of physical accuracy, it allows for all the most 
important physical ingredients that may carry out a role in 
the evolution of a spherical cluster. These include, among 
others, self-gravity, two-body relaxation, stellar evolution, 
collisions, binary stars etc and, undoubtedly, the interaction 
with a central BH and the role of a mass spectrum. The 
specific advantage of the so-called "gaseous model" to other 
simulations methods (to be presented in 2.2) is that the 
simulations are comparatively much faster , since they are 
grounded on numerical integration of a relatively small set of 
partial differential equations with just one spatial variable, 
the radius r. In addition, all quantities of interest are acces- 
sible as smooth functions of r and time and this allows one 
to investigate in detail clean-cut aspects of the dynamics 
without being hindered by the important numerical noise 
particle-based methods (TV-body and Monte Carlo) suffer 
from. 

In this paper we concentrate on the simplest version of 
the gaseous model which includes the interaction between 
a central BH and its host cluster. In particular, we assume 
that all stars are sun-like, neglect stellar evolution, direct 
collisions between stars and the role of binary stars. Also, 
the only interaction between BH and the stellar system is 
tidal disruptions (besides the BH's contribution to the grav- 
itational field) , and we undertake that the BH stays fixed at 
the centre of the cluster. While admittedly very simplified, 
we reckon this idealised situation warrants consideration. 
First, it helps us to establish that the gaseous model is also 
able to treat more complex situations in phase space than 
self-gravitating star clusters, such as those caused by loss- 
cone accretion and a central BH; second, experience with the 
gaseous model and other methods showed us how intricate 
the interplay between various physical effects can become 
during the evolution of clusters and so we feel compelled to 
first consider the simplest models to develop a robust un- 
derstanding of the mechanisms at play. 

The structure of this paper is as follows: In the sub- 
section 2.1, we explain the physics of the problem. In 2.2 
we introduce briefly the various analytical and numerical 
methods used so far to investigate spherical clusters with 
central BH and summarise their key results. Later, in sec- 
tion 3, we present the gaseous model and, in particular, how 
the interaction between the stars and the BH is included. 
The interested reader may find in the Appendix A a more 
technical description of the numerical method used to solve 
the equations of the gaseous model. Section 4 is devoted to 
results obtained in a first set of simulations. Finally, in 5, 
we draw conclusions about this first use of this new code 
version and present what our future work with it will likely 
consist of. 



2 LOSS-CONE ACCRETION ON TO MASSIVE 
BHS: THE DIFFUSION MODEL 

2.1 Previous theoretical and numerical studies 

If we arrange numerical methods for stellar dynamics in 
order of validity and both increasing the spatial resolu- 
tion and decreasing the required computational time, we 
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can distinguish four general classes. Th e most direct a p- 
proach is the so-called TV-body method jAarsethlll999allbT: 
ISpurzemlll99flr i . Monte Carlo codes are also particle-based, 
but rely on the assumptions that the system is spheri- 
cally symmetric and in dynamical equilibrium and treat 
the relaxation in the F okker-Planck approx i mation (see 3.2 



the relaxation in the ^okkcr-Flanck approximation (sec 6.2) 
jFreitag fc Benzll2004l; iFregeau et alj|2003l: iFreitag fc Beiiz 



20011: iGiersz fc Spurzeml 120031 : iGiersd Il998t IJoshi et al 
20011) . Ensuingly, we have the two-dimensional numerical 



direc t solut i ons of the Fokker-Planck equation (iTakahashil 
119971. Il99fj. Il995l) . and the gaseous model s. The idea 
of this model goes back to lHachisu et al.l lll978l) and 
iLvnden-Bell fc Eggletonl il980f) . who first proposed to treat 
the two-body relaxation as a transport process like in a 
conducting plasm a. They had been developed fu r ther by 
IBettwieserl Il983tl: IfBettwieser fc Sugimotol (Il984h : iHeggid 
I^^^ T lHe^ie^ r^miMn6mirir i989h. T heir present form, 
publis h ed in iLou is^^P urzerrJ lll99lT) : IGiersz fc Spurzeml 
ill994f ); rSpurzem^^Tak^iashirill995l) improves the detailed 
form of the conductivities in order to yield high accuracy (for 
comparison with TV-body) and corre ct multi-mass m odels. 
Thi s point has bee n made already in|Srjurzem| jll992l) . 

IPeeblesj <1972T): IShapiro fc Lightmanl il97rJl and e spe- 
cially iFrank fc Reesl TT976D and iBahcall fc Wold (Il976l) ad- 
dressed the problem of a stationary stellar density profile 
around a massive star accreting BH. They found that, un- 
der certain conditions, the density profile p oc r~ 7//4 is estab- 
lished in the region where the BH's gravitational potential 
well dominates the self-gravity of the stars 2 . 

The problem of a star cluster with a massive cen- 
tral star-accreting BH has been widely coped with Fokker- 
Planck numerical models. This approach was useful in or- 
der to test the solidness of the method to reproduce the 
p oc r~ 7 / 4 stationary density profil e, since loss-cone ac- 
cretio n disturbs such a density cusp dOzernoi fc Reinhardtl 
I1978T) . The authors show that the stationary density profile 
follows from their stellar-dynamical equation of heat trans- 
fe r by scaling arguments which are analogous to those given 
in lShaoiro fc Lightmanl (Il97rj) . 



2.2 Loss-cone physics 

We can express the tidal radius in terms of the 
internal stellar structure re-w riting equation (2) of 
lAmaro-Seoane fc Spurzeml <l200lT l (in pursuit AS01) 



r t oc 



\ 1/3 

■Mbh ) 

X b TTp I 



(1) 



where Mbh denotes the mass of the central BH, p the mean 
stellar internal density, n is the polytropic index (stars are 
supposed to be polytropes) and Xb is a parameter propor- 
tional to the gravitational binding energy of the star that 
describes effects of the internal stellar structure. We assume 
that a star is disrupted by tidal forces when it crosses the 
tidal radius. The free parameter e a a (accretion efficiency) 



determines the mass fraction of the gaseous debris being ac- 
creted on to the central BH (e c ft = 1 corresponds to 100% 
efficiency) . 

There are two concurrent processes driving stars to- 
wards the tidal radius; namely the energy diffusion and the 
loss-cone accretion. In the first case, stars on nearly circular 
orbits lose energy by distant gravitational encounters with 
other stars and in the process their orbits get closer and 
closer to the central BH. The associated energy diffusion 
time-scale can be identified with the local stellar -dynamical 
relaxa tion time generalised for anisotropy as in IBettwieserl 
Jl983ft: 



9 



.( fft 2 /2) 



160FG 2 m*p*(r)lnA' 



(2) 



Here oy, at are the radial and tangential velocity dispersions 
(in case of isotropy 2o> = a 2 ), p*(r) is the mean stellar mass 
density, N the total particle number, G the gravitational 
constant, m* the individual stellar mass and 



In A = In (p m ax/po) = In ( 7 AT) (3) 

is the Coulomb logarithm. We set 7 = 0.11 iGiersz fc Heggiel 
Il994l) . In this expression p max is an upper limit of p, the 
impact parameter, po is the value of p that corresponds to 
an encounter of angle ip — 7r/4, where tp — (tt — £)/ 2 if e: is 
defin ed to be the deflection angle of the encounter iSpitzerl 
Il987l) . In the vicinity of the BH (r < rh, see Sec . l3~5)l . one 
should be aware tha t A sa Mbh/m* jBahcall fc Woljll976l : 
iLightman fc Shapiro! Il977ll but, for simplification, here we 
shall use Eq.|S| (strictely speaking only valid at distances 
r > rh) everywhere. 

For a more detailed discussion of the energy diffusion 
process a nd its description in the cont ext of the moment 
model see iBettwieser fc Spurzeml dl986l) . 

As regards the second process, the loss-cone accretion, 
stars moving on radially elongated orbits are destroyed by 
tidal forces when they enter the tidal radius r t . A star will 
belong to the loss-cone when its peribarathron 3 (distance of 
closest approach to the BH, see Fig0 is less than or equal to 
the tidal radius r t , provided that its orbit is not disturbed 
by encounters. Thus, the loss-cone can be defined as that 
part of stellar velocity space at radius r, which is given by 



\v t \ < vi c (r) 



rt 



• ^2[4>{r t )-<t ) {r)]+v r {rY (4) 



(see AS01). In the last formula v r , v t are the radial and 
tangential velocity of a star and 

4>(n) *(r) = ^ + Mn) - ^ - Mr) (5) 

At distances r > r t we can approximate this expression 
taking into account that GMbh/n 4>{r) and (f>(rt), 



2 We must mention h e re the legwork done twelve years before this 
analysis by iGurevichl Il964f) . since he got an analogous solution 
for the distribution of electrons in the vicinity of a positively 
charged Coulomb centre. 



3 This word fits quite well the idea of closest approach to this 
"sinking" hole for it has the meaning of no return. The barathron 
(P&pctOpov) was in the ancient Greece a cliff down to an un- 
reacheable or unseen place where criminals were thrown. 



4 P. Amaro-Seoane, M. Freitag and R. Spurzem 



Vtg (max) 




the three-dimensional velocity space is re-populated within 
a time- scale 



where 



o 



tiii — tj- 



fd\/ 



fd 3 v; 



(9) 



(10) 



the subscript oo denotes an integration over the velocity 
space as a whole, and the subscript Ic means that the inte- 
gration over the loss-cone part of the velocity space is given 
by Eq. 0. As a matter of fact, Q can be envisaged as the 
fraction of the surface of a velocity ellipsoid which is cut 
out by the loss-cone. Close to the tidal radius r t , and for 
appreciable amounts of stellar-dynamical velocity dispersion 
anisotropics, our method desc ribes the loss-cone size Q more 
exact ly than those given by (iFrank fc Reeslll976l : Ida Costal 
Il98lf) . On the other hand, their models can be "recovered" 
in the limit of r > r t and isotropy, 2a 2 . ~ a 2 (denoted as 
the "small loss-cone approximation"), where 



Figure 1. Definition of the peribarathron as the distance of closest 
approximation of the star in its orbit to the BH. In this point 
the radial component of the velocity of the star cancels and the 
tangential component is maximum. In the figure "r p " stands for 
the peribarathron radius and "rt" for the tidal radius 



This means that 



vi c (r) 



r t 2GMy 



n. 



(6) 



On the other hand, 



a r (r) = 



GM bh 



Thus, it is in fair approximation 
vic{r) 



-1/2 



(7) 



(8) 



For a deeper a nalysis on l oss-co ne phenomena see AS01. 

Similar to Ida Costal l|l98ll ). we define two time-scales: 
tout to account for the depletion of the loss-cone and tin for 
its replenishment. For a BH, t ou t is equivalent to one crossing 
time, since it is assumed that a star is destroyed by just one 
crossing of the tidal radius. In AS01 we mentioned that in 
the special case that the central object is a super-massive 
star, t ou t = nt cross . Here n > 1 is the number of passages 
until a star is trapped in the central object. 

The loss-cone is replenished by distant gravitational en- 
counters that change the angular momentum vectors of the 
stars. To estimate t ln , we make the assumption that ran- 
dom gravitational encounters thermalise the whole veloc- 
ity space at some given radius r after a time-scale of the 
order of t re iax- As a first approximation, the fraction SI of 



n 



2 / 2 _ / 



(11) 



This is equivalent to their definition of 9\ c if Q = 9 2 c /4l is 
adopted. 

Where the loss-cone effects can be neglected, a Schwarz- 
schild-Boltzmann type distribution function can be as- 
sumed, 



/ oc exp 



(V T - {v r }) 2 v\ 



2al 



(12) 



Third order moments of the velocity distribution that rep- 
resent the stellar-dynamical energy flux do not alter such a 
distri bution function significantly feettwieser fc Sugimotol 
119851) . 

An important quantity is the critical radius r cr it. Let 
#o be the average quadratic deflection angle produced by 
relaxation during t out (= tcross here), 



Q 2 tout 



trclax 



Then, by definition, 



?lc(r C rit) = fc(r C rit), 



(13) 



(14) 



which is equivalent to tout(fcrit) = 4ti (r C rit)- For most clus- 
ters where such a radius can be defined, 0\ c > On inside r cr i t 
while the opposite holds outside. This means that, at large 
radii, relaxation is efficient enough to make stars diffuse into 
and out of loss cone orbits over a time scale t ou t so that the 
distribution function is not appreciably depleted in the loss 
cone. Conversely, deep inside r cr it, the loss cone orbits are 
essentially empty and the flux of stars into this domain of 
phase-space (and into the BH) can be treated as a diffusive 
process because the size of one individual step of the ve- 
locity random walk process, <9d, is (much) smaller that the 
characteristic size of the problem, 0i c . 

Note that a critical radius does not necessarily exist (see 
AS01). For instance, if one assumes that gravity of the BH 
dominates the stellar self-gravitation and that the density 
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profile follows a power-law, p oc r~ a , one has 8f c oc r -1 , 
#o oc r 3 ~ a and a critical radius would not exist for x > 4. 

Now we want to generalise the stationary model (see 
AS01), which assumes an empty loss-cone within r cr it and 
a full loss-cone elsewhere, by means of a simple "diffusion" 
model, which is derived from the above considerations; this 
means that the filling degree of the loss-cone K can be con- 
tinuously estimated within its limiting values, 



K G [0, 1] 



(15) 



Let / be the unperturbed velocity distribution (without loss- 
cone accretion); if the loss-cone is empty and we neglect 
the angular momentum diffusion, / = inside the loss-cone 
(and unchanged elsewhere in the velocity space). In point of 
fact, / will have a continuous transition from nearly unper- 
turbed values at large angular momenta to a partially de- 
pleted value within the loss-cone. This value is determined 
by the ratio of tin and t ou t . Such a smooth transition of the 
distribution function given as a function of angular momen- 
tum f(J) has been derived from self-consistent models of 
angular momentum diffusi on (e.g. ICohn fc Kulsrudll 19781 or 
iMarchant fc Shapirdll98(I) . We approximate f(J) by a dis- 
tribution function that has a sudden jump just at the value 
Jmin = m*v\ c from an unperturbed value fo given by the mo- 
ment equations (assuming a Schwarzschild-Boltzmann dis- 
tribution) to the constant lowered value 



/ = Kf , with < K 1 



(16) 



within the loss-cone (i.e. J < J m i n or |i>t| < v\ c ). This implies 
that, as means to compute the mean mass density of loss- 
cone stars, we have to calculate the integral 



Pic 



Kf d%. 



(17) 



As we assume relaxation is due to a large number of 
small-angle deflections and can thus be seen as a diffusive 
process in velocity space, the probability P(6) that a star is 
scattered in an angle 9 in a time t out is 



P{8) = 



cxp(-0 2 /^), 



The distribution is normalised to one, 



P{6)d6 : 



(21) 



(22) 



and has the property that its mean square value is d 2 D . A 
star remains in the loss-cone during a time t ou t if its RMS 
diffusion angle is smaller than the loss-cone angle 0i c . The 
probability for this to happen is 



Plc = 



P{9) d6 = erf (^4 t in /to 



(23) 



In the case that a star is unperturbed by the rest of the 
stellar system, it will sink on to the central BH in a time 
tout- Actually, this is a somehow simplified description of the 
physical process, for part of the loss-cone stars will be scat- 
tered out of it before they slump. The required time for this 
event is tin, since in this time-scale the angular momentum 
vector will change (due to distant encounters) on an amount 
that is comparable with the size of the loss-cone in the an- 
gular momentum space. For this reason we have introduced 
the quantity Pi c in Eq. 12UH . 

Bluntly speaking, the effective time-scale that describes 
the loss-cone depletion allowing for perturbation due to an- 
gular momentum diffusion is 



And then, accordant with the definition of 51, 



tout, eff — tout I ' P\c 



(24) 



Pic, full = pfi, 



(18) 



in the case that we have a full loss-cone. 
In regard to the radial and tangential stellar velocity disper- 
sions in the loss-cone cric,r and <ti c , t, we can compute them 
using second moments integrated over the loss-cone part of 
velocity space. As for the definition of the quantities E r and 
E t used in Sect. 3, 



2 

°"lc, r 
2 

0"lc,t 



= E r a', 
= E t er? 



(19) 



in the small loss-cone approximation we have that E r ~ 1 
and E t < 1. 

The arguments about the time-scales that have led us 
to the derivation of ti n and t ou t guide us also to the follow- 
ing diffusion equation for the time evolution of the spatial 
density pi c = KpQ, of loss-cone stars: 



rfpic 
dt 



Pic Pic P^ 

^out 



Pic 



tin 



(20) 



In this equation, the second term on the right hand is the 
refilling term due to relaxation. 



As a matter of fact, this definition ensures us that far outside 
of the critical radius the loss-cone depletes in a time that 
grows infinitely, as it is physically expected. In the regime 
where r <C r cr it , Pic tends asymptotically to 1 and to where 
r ^> r cr it, passing through a transition zone at r = r cr it- 

We can consider Eq. (1201 as an ordinary differential 
equation for K = pi c /(pfi) if we assume that the stellar 
density and the loss-cone size are time-independent. Trans- 
port phenomena can be neglected, for they are related to the 
relaxation time, and f ro i ax 3> tin, tout- Bearing this in mind 
we can get an analytical solution K(t) for the differential 
equation with the initial condition that K(t)\ t(l = Kq, 



^(t) = A-oexp(- Pl ^, ( ^ to) )4 



tout 



tout 
Pic tin £ 



1 — exp 



Plc£(t-t ) 



(25) 



In the last equation we have defined £ for legibility reasons 
as follows, 



f := 1 + (tout /Pic tin). 



(26) 



6 P. Amaro-Seoane, M. Freitag and R. Spurzem 



For r = r cr it, with t- m = tout, the stationary filling degree of 
the loss-cone turns out to be 

K x := lim K(t) = J. (27) 

t — »oo z 

Note that lMilosavlievic fc Merrittl (I2003T) recently gave a de- 
tailed summary of loss-cone effects. They derived expres- 
sions for non-equilibrium configurations. They employ a 
rather different treatment for the diffusion since they tackle 
the problem of binary BHs scattering. 



3 THE GASEOUS MODEL 
3.1 Introduction 

In this section we will introduce the fundamentals of the 
numerical method we use to model our system. We give a 
brief description of the mathematical basis of it and the 
physical idea behind it. The system is treated as a contin- 
uum, which is only adequate for a large number of stars and 
in well populated regions of the phase space. We consider 
here spherical symmetry and single-mass stars. We handle 
relaxation in the Fokker-Planck approximation, i.e. like a 
diffusive process determined by local conditions. We make 
also use of the hydrodynamical approximation; that is to 
say, only local moments of the velocity dispersion are con- 
sidered, not the full orbital structure. In particular, the ef- 
fect of the two-body relaxation can be modelled by a local 
heat flux equation with an appropriately tailored conductiv- 
ity. Neither binaries nor stellar evolution are included at the 
present work. As for the hypothesis concerning the BH, see 
3.4. 



3.2 The local approximation 

We treat relaxation like the addition of a big non-correlated 
number of two-body encounters. Close encounters are rare 
and thus we admit that each encounter produces a very small 
deflection angle. Thence, relaxation can be regarded as a 
diffusion process 4 . 

A typical two-body encounter in a large stellar system 
takes place in a volume whose linear dimensions are small 
compared to other typical radii of the system (total system 
dimension, or scaling radii of changes in density or velocity 
dispersion). Consequently, it is assumed that an encounter 
only changes the velocity, not the position of a particle. 
Thenceforth, encounters do not produce any changes Ax, 
so all related terms in the Fokker-Planck equation vanish. 
However, the local approximation goes even further and as- 
sumes that the entire cumulative effect of all encounters on a 
test particle can approximately be calculated as if the parti- 
cle were surrounded by a very big homogeneous system with 
the local distribution function (density, velocity dispersions) 
everywhere. We are left with a Fokker-Planck equation con- 
taining only derivatives with respect to the velocity vari- 

4 Anyhow, it has been argued that rare deflections with a 
large angle may play a important role in the vicinity of a BH 
jLin fc TremainJl98dl . 



ables, but still depending on the spatial coordinates (a local 
Fokker-Planck equation). 



3.3 A numerical anisotropic model 

For our description we use polar coordinates, r 9, <f>. The 
vector v = — r,9,4> denotes the velocity in a local 

Cartesian coordinate system at the spatial point r, 6, 4>- For 
succinctness, we shall employ the notation u = v r , v = Vg, 
w = v$. The distribution function /, is a function of r, t, u, 
v 2 + w 2 only due to spherical symmetry, and is normalised 
according to 



p(r,t) = / f(r,u,v 2 +w 2 ,t)dudvdw. (28) 



Here p(r, t) is the mass density; if m* denotes the stellar 
mass, we get the particle density n = p/m*. The Euler- 
Lagrange equations of motion corresponding to the La- 
grange function 



C = - (r 2 + r 2 9 2 + r 2 sm 2 8 2 ) - $(r, t) (29) 
are the following 



u — — 

v — — 
w = — 



9$ 


+ 


2 i 2 
V +W 


dr 


r 


uv 
r 


+ 


w 2 
rtanO 


WW 




vw 


r 




rtanO 



(30) 



And so we get a complete local Fokker-Planck equation, 



df df . df . df . df ( 5 A 

i +v ^ +Vr oi +ve d!r g +v ^ = {i) Fp (31) 

In our model we do not solve the equation directly; we 
use a so-called momenta process. The momenta of the veloc- 
ity distribution function / are defined as follows 



/■ + oo 

<i,j,k>:= / v 1 r viv^f(r,v r ,ve,v < f > ,t)dVrdvodv ct y; (32) 



We define now the following moments of the velocity distri- 
bution function, 
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< 0,0,0 >:= 

< 1,0,0 >:= 

< 2,0,0 >:= 

< 0,2,0 >:= 

< 0,0,2 >:= 

< 3,0,0 >:= 

< 1,2,0 >:= 

< 1,0,2 >:= 



p — J fdudvdw 
u — u fdudvdw 



p r + pu 2 — / u 2 fdudvdw 



pe — v fdudvdw 



(33) 



P4> = w fdudvdw 



F r + 3up r + it 3 = / v? fdudvdw 



Fe + upe = J uv 2 fdudvdw 
Ftf, + up,f, = j uw 2 fdudvdw, 



where p is the density of stars, u is the bulk velocity, v r 
and vt are the radial and tangential flux velocities, p r and 
pt are the radial and tangential pressure s, F r is the radial 
and F t the tangential kinetic energy flux l|Louis fc Spurzeml 
Il99ll) . Note that the definitions of pi and Fi are such that 
they are proportional to the random motion of the stars. 
Due to spherical symmetry, we have pg = p$ =: pt and 
Fe = Ftf, —: Ft/2. By p r = pa 2 and pt = pa 2 the random 
velocity dispersions are given, which are closely related to 
observable properties in stellar clusters. 
F — (F r + Ft)/ 2 is a radial flux of random kinetic energy. In 
the notion of gas dynamics it is just an energy flux. Whereas 
for the 9— and <f>— components in the set of Eqs. 1331 are 
equal in spherical symmetry, for the r and t- quantities this 
is not true. In stellar clusters the relaxation time is larger 
than the dynamical time and so any possible difference be- 
tween p r and p t may survive many dynamical times. We 
shall denote such differences anisotropy. Let us define the 
following velocities of energy transport: 



v t = 



3p r 

Ft_ 

2p t 



+ u, 
+ u. 



(34) 



In case of weak isotropy (p r ~Pt) 2F r = 3F t , and thus v r — v t , 
i.e. the (radial) transport velocities of radial and tangential 
random kinetic energy are equal. 

The Fokker-Planck equation 131H is multiplicated with 
various powers of the velocity components u, v, w. We get 
so up to second order a set of moment equations: A mass 
equation, a continuity equation, an Euler equation (force) 
and radial and tangential energy equations. The system of 
equations is closed by a phenomenological heat flux equation 
for the flux of radial and tangential RMS (root mean square) 
kinetic energy, both in radial direction. The concept is phys- 
ically similar to that of lLvnden-Bell fe Eggletonl il980f> . The 
set of equations is 



dp . 1 d , 2 s n 
at r 2 dr 



du 
dt 



GM r 1 dpr 
r 2 p dr 



Pr - pt 
pr 



= 



dp r ld.o \ „ du 19,2,-, 

—Li. _| (r up r ) + 2p r | (r F r 

dt r 2 dr T T dr r 2 dr r 



(35) 



2F t 



4 (2 Pr - Pt) 

5 AA^rolax 



dpt 19,2 s . „ Pt u 1 d . 2 „ > 
-dt+^d-r {r UPt)+2 — + 2^a-r {rFt) 
+ Ft = 2 (2p r -pt) 

r 5 A^rclax 

where Aa is a numerical constant related to the time-scale 
of collisional anisotropy decay. The value chosen for it has 
been discussed in comparison wi th direct simulations per - 
formed with the iV-body code jGiersz fc Spurzeml ll994Tl . 
The authors find that Aa = 0.1 is the physically realistic 
value inside the half-mass radius for all cases of N, provided 
that close encounters and binary activity do not carry out 
an important role in the system, what is, on the other hand, 
inherent to systems with a big number of particles, as this 
is. 

With the definition of the mass M r contained in a sphere of 
radius r 



dM r 
dr 



= 47rr 



(36) 



the set of Eqs. il-'-iot is equivalent to gas-dynamical equations 
coupled with the equation of Poisson. To close it we need an 
independent relation, for moment equations of order n con- 
tain moments of order n + 1. For this intent we use the heat 
conduction closure, a phenomenological approach obtained 
in an analog ous way to gas dynamics . It w as used for the 
first time bv iLvnden-Bell fc Eggletonl lll980F) but restricted 
to isotropy. In this approximation one assumes that heat 
transport is proportional to the temperature gradient, 



F = -k— = -A— 
dr dr 



(37) 



That is the reason why such models are usually also called 
conducting gas sphere models. 

It has been argued that for the classical approach A oc 
A 2 /t, one has to choose the Jeans' length Aj = a 2 /(4ivGp) 
and the sta ndard Chandrasekhar local r elaxation time 
irciax oc cr 3 /p iLvnden-Bell fc Eggletonll980Tl . where A is the 
mean free path and r the collisional time. In this context we 
obtain a conductivity A oc p/o. We shall consider this as 
a working hypothesis. For the anisotropic model we use a 
mean velocity dispersion a 2 — (o 2 + 2a 2 )/ 3 for the temper- 
ature gradient and assume v r = vt iBettwieser fc Spurzeml 
Il986f) . Forasmuch as, the equations we need to close our 
model are 



X 



da 2 



47rGpt rc i ax dr 



= 



(38) 
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3.4 Inclusion of the central BH in the system 



In this subsection we discuss the way we cope with the loss- 
cone in our approach. For this aim we accept the following: 

1. The system has central (r — 0) fixed BH 

2. Stars are totally destroyed when they enter r t 

3. Gas is completely and immediately accreted on to the BH 

As regards the first point, one should mention that the 
role of brownian motion of the central BH can be important; 
as a matter of fact, for a cluster with core radius i? C ore, 
cquipartition predicts a wandering radius R wa , n of order 



R* 



Re 



y/M<./Mbh, 



(39) 



which is larger than the tidal disruption radius for 
BHs less massive than 10 9 M W if the core radius 
is 1 pc teahcall fc Woll 1 19761: iLin fc Tremaind Il980t 
IChatterle^et^dT2002T l~The wandering of a MBH at the cen- 
tre of a cuspy cluster has been simulated bv lDorband et all 
l|2003h with a JV-body code allowing N = 10 6 . They find 
that RMS velocity of the quickly reaches equipartition with 
the stars but do not comment on the wandering radius. In 
Appendix |0 we present a simple estimate suggesting that, 
in a cusp p oc r~ a , the wandering radius may be much re- 
duced, Ewan oc a(m*/A4bh) 1/ '' 2_c *', where a is typical length 
scale for the central parts of the cluster. For a ^ 1.5, one 
would then expect 7? wan to be smaller than R t for black 
holes as light as 2000 Mq , but this arguments neglects the 
flattening of the density profile due to loss-cone accretion. 
Further iV-body simulations are clearly required to settle 
the question and, in particular, if i? wan is larger than R t , 
to establish the effect of the motion of the MBH on dis- 
ruption rates, which can be either increased or decreased 
jMagorrian fc Tremainelll999l) . 

We arrogate that at any sphere of radius r the transport 
of loss-cone stars in the time-scale £ ou t, e ff towards the centre 
happens instantaneously compared with the time step used 
for the time evolution. 

Hence, the local density "loss" at r is 



<5p\ Pic Pic , 

St ) , tout 

with pic = Kttp. 

This corresponds to a local energy loss 5 of 



(40) 



8p erf. 

Spa? 
St 



Sp\ „ 2 



(41) 



E r and Et are worked out integrating over the velocity dis- 
tribution part that corresponds to the loss-cone with the 
approximation u <C ov and v\ c <JC at, 



5 By loss we mean here transport of mass and kinetic energy 
toward the central BH, for it is lost for the stellar system. 



E r 

E, 



V\c/<7t < 1 



(42) 



Thereupon, the mass accretion rate of the central BH can 
be calculated as 



M = 



-eefl 



Rto 



^ | 4irr 2 dr. 

<W,c 



(43) 



Here 7? tot stands for the total radius of the stellar system. 
The accretion efficiency has been set throughout our calcu- 
lati ons to e e ff = 1; for a discu ssion on different e e g -values 
see lMarchant fc Shapiro! Jl98tf) . 

The complete set of Eqs. 1)3 5 ^ including the local ac- 
cretion terms of the type (S/8t)\ c for energy diffusion and 
loss-cone accretion are solved implicitly. For every time step 
the mass of the BH and the filling degree K of the loss-cone 
are brought up to the new state of the system. The time 
step is chosen in order to keep the maximum changes of the 
variables below 5% . 

For the model calculation we have utilised for the boundary 
conditions that at the outer limit, Rtot = 10 4 pc, we impose 
u = F — and M r = M to t- No stellar evaporation is allowed. 
At the centre, the usual boundary conditions for the gaseous 
model are u — M T — F = but the central point ri = is 
not explicitly used when there is a BH, for obvious reasons. 
Instead, one imposes that all quantities vary as power-laws, 
d In x/d In r — C Bt inside the first non-zero radius of the 
discretisation mesh, T2 = 1.7 x 10 -6 pc (see Appendix A). 

3.5 Units and useful quantities 

The units used in the computations correspond to the so- 
called iV-body unit system, in which G = 1, the total ini- 
tial mass of th e stellar clus t er is 1 and its initial tot al en- 
ergy is -1/2 <Henodll97ll iHeggie fc Mathieul Il986tl . For 
the simulations presented here, the initial cluster structure 
corresponds to the Plummer model whose density profile 
is p(r) = po(l + (r/Rpi) 2 ) 5 ^ 2 , where _Rpi is the Plummer 
scaling length. For such a model the iV-body length unit is 
Ui = 16/(3tt) R Ph 

In the situations considered here, the evolution of the 
cluster is driven by 2-body relaxation. Therefore, a natural 
time scale is the (initial) half- mass relaxation time. We use 
the definition of ISpitzerl (I19871) . 



T rh (0) 



0.1387V 
In A 



Rl 



/2 



GMc 



1/2 



(44) 



For a Plummer model, the half-mass radius is R1/2 = 
0.769Wi = 1.305 Rpi. A4 c i is the total stellar mass. 

For a cluster containing a central BH, an important 
quantity is the influence radius, enclosing the central region 
inside of which the gravitational influence of the BH dom- 
inates over the self-gravity of the stellar cluster. The usual 
definition is r n = G7Mbh/o"o> where 00 is the velocity dis- 
persion in the cluster at a large distance from the BH. As 
the latter quantity is only well defined for a cluster with 
an extended core, we use here the alternate and approxi- 
mate definition M r (r n ) = A^bh, he. r^ is the radius of that 
encloses a total stellar mass equal to the mass of the BH. 
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4 RESULTS 

We study the evolution of a stellar cluster with a so-called 
"seed BH" at its centre. We consider two possible con- 
figurations for the stellar system; one of a total mass of 
M tot = 10 5 M Q and another of 10 6 M o . For the initial BH 
mass, we have chosen .Mbh(O) = 50Mq and 500Mq and we 
model it as a Plummer of Rpi — lpc. Even though it would 
be more realistic to set the efficiency parameter e c g = 1/2, 
we choose here e c a — 1 for historical reasons. Nonetheless, 
here we study additionally the influence of a variation of the 
stellar structure parameter Xb, since it influences the tidal 
radius and hence the accretion rates (see Eq.0. For this in- 
tent we compare case a (xb = 1) with another one in which 
we choose the value Xb = 2 (case b). As regards the physi- 
cal meaning, the stars of case b have twice as much internal 
binding energy than case a. 

The cluster evolves during its pre-collapse phase up to a 
maximum central density from which the energy input due 
to star accretion near the tidal radius becomes sufficient 
in order to halt and reverse the core collapse. Immediately 
afterwards, the post-collapse evolution starts. At the begin- 
ning of the re-expansion phase, the BH significantly grows to 
several 10 solar masses. Thereon, a slow further expansion 
and growth of the BH follow. 

In Fig.|21 we follow the evolution of the mass of a central 
BH in a globular cluster of 10 stars of 1 Mq . Panel (a) 
shows the mass of the BH as function of time. On panel (b) , 
we present the accretion rate on to the BH, i.e. its growth 
rate. For jMbh(O) = 50 Mq, the early cluster's evolution is 
unaffected by the presence of the BH which starts growing 
suddenly at the moment of deep core collapse, around T ~ 
14.5 T r h(0). In Fig.[3]we follow the same evolution for the 
case of a stellar cluster of 10 6 stars. 

From Figs.|2] and |3| we can see that the differences be- 
tween the cases a, b and c are nearly negligible after core 
collapse. In general, the structure of the cluster at late times 
is nearly independent of A / fbh(0) and x\,. From these plots 
we can infer that this occurs since core collapse leads to 
higher densities if the initial BH mass is smaller and thus 
the integrated accreted stellar mass increases. 

We exhibit the evolution of the structure of the clus- 
ter for case a with 10 s stars in Fig.|3] With dotted lines we 
plot various Lagrangian radii, for mass fractions ranging be- 
tween 10 -3 % (which formally corresponds to only one star) 
to 90 %. Only the mass still in the stellar component at a 
given time is taken into account. Moreover, the evolution of 
the influence radius (solid line, defined as the radius enclos- 
ing a stellar mass equal to the BH mass) and critical radius 
(dashed line) are shown, so that one can infer the percentage 
of the stellar mass embodied within them at a certain mo- 
ment. For late time, one obtains self-similar evolution with 
size increasing like R oc T 2 ^ 3 , as expected for a system in 
which the central object has a small mass and the energy 
production is confined to a small centra l volume. (iHeno nl 
Il965l: IShapirolll977l: iMcMillan et~aT1ll98ll lGoodmanlll984l) . 
We consider too the case of a 10 6 stars in Fig.|S| for which 
the R oc y 2 / 3 expansion is a poor approximation because, 
at late times, the BH comprises of order 40 % of the system 
mass. 

We observe in Fig.|S|that for the evolved post-collapse 
model the spatial profile of the stellar density has a power 
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Figure 2. Evolution of the mass of a central BH in a globular 
cluster of 10 5 stars of 1 Mq. We considered three cases. In case a 
(solid line), the initial BH mass is .Mbh(O) = 50 Mq and Xb = 1, 
case b (dashes) has the same initial BH mass but Xb = 2 while 
case c (dash-dot) corresponds to jVfbh(O) = 500 Mq and Xb = 1. 
An accretion efficiency of e c ff = 1 is assumed. Panel (a) shows the 
mass of the BH as function of time and panel (b) the accretion rate 
on to the BH. At late times, the mass of the central BH increases 
like A4bh °c T — 1,2 as predicted by simple scaling arguments (see 
text). 



law slope of p oc r-~ 7//4 in the region r C rit < r < rh, where 
r\ is the influence radius. The density profile flattens for 
t < r cr i t due to the effective loss-cone accretion. For the 
same post-collapse moment we display the surface density 
for case a in Fig.E] 

iLightman fc Shapirol (119771) proved that within r C rit , accord- 
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Figure 3. Same as Fig.|2]but for a cluster of 10 6 stars with the 
same size (Rp = lpc). 



ing to their Eq. (71) (where they assume Q -C 1, small loss- 
cone), q will continously vary from 1.75 to -co. 

Regarding the three dimensional velocity dispersion, in 
Fig.Qwe can see coming a slope of -1/2 at the inner region; 
but extending below the one-star radius does not make much 
sense. A slope of -1/2 is what one would expect from Ke- 
pler's third law and a simple application of Jeans equation, 
with the assumptions that (1) dynamical equilibrium holds, 
(2) the gravity is dominated by the central BH, (3) the den- 
sity follows a pure power-law and (4) the anisotropy a t /a r , 
is constant, indeed predicts a cx r^ 1 ^ 2 . In Appendix |E| 
we show that the Jeans equation for stationary equilibrium 
actually describes the central regions of the cluster quite 
well. The reason why the velocity dispersion does not follow 
closely the "Keplerian" profile has to do with the fact that 
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Figure 4. Evolution of the radii of spheres enclosing the indicated 
fraction of the total mass, Lagrangian radii, for case a. The mass 
fractions range from 10 -3 % to 90 %. The influence and critical 
radius arc displayed (solid and broken line). See text for further 
explanation. 
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Figure 5. Same as Fig.lll but for a cluster of 10 6 stars. 



none of assumptions (2)- (4) exactly holds all the way from 
the influence radius inward. 

Figures lTITl and IH1 give the plots of the projected density 
and velocity dispersions for the late post-collapse model. 

The effects of anisotropy are studied in Fig.|H| where we 
can see that the external parts of the cluster are dominated 
by radial orbits. Inside the critical radius (indicated by a 
star symbol), one notices a slight tangential anisotropy, an 
effect of the depletion of loss-cone orbits. At large radii, the 
velocity distribution tends to isotropy as an effect of the 
outer bounding condition imposed at 10 4 pc. 
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Figure 6. Density profile for a 10 Mq globular cluster with 
•Mbh(O) = 50 Mq (case a) at 10 Gyrs. The round dot indicates 
the influence radius (M r = A^bh)> the star the critical radius and 
the square the radius below which the description of the cluster 
as a continuum loses significance because the enclosed mass is 
smaller than lMg. As expected, for the zone between the crit- 
ical and influence radii, the density profile closely reassembles 
a power-law of exponent —7/4. At that stage, the structure of 
case c (A4bh(0) = 500 Mq) is extremely similar. We can see that 
from the "1-star" radius onwards the slope of the curve shows a 
tendence to increase. This is due to the fact that the loss-cone 
size is artificially limited for stability purposes. On the other 
hand, the slope comprised between the critical radius and the 
1- star radius is consistent wit h the arguments given in the work 
of lLight man fc Shapiro! <1977l) (see text for further explanation). 



The effects of anisotropy in the stellar system can also be 
seen in Fig. llQI where we plot the components along the 
line of sight ctlos (solid line) and on the sky ("proper mo- 
tions"). The latter is decomposed into the radial direction 
(i.e. towards/away from the position of the cluster's centre) 
component, apM,r (dashes) and the tangential component, 
o"PM,t (dash-dot). Note that the radial anisotropy in the out- 
skirts of the cluster reveals itself as a radial "proper motion" 
dispersion slightly larger than the other components. For 
an isotropic velocity dispersion, all three components would 
be equal. Despite loss-cone effects, there is no measurable 
anisotropy at the centre. 

The loss-cone induced anisotropy could be detected only 
if one could select the stars that are known to be spatially 
close to the centre (and not only in projection) as would be 
feasible these stars happened to be of a particular popula- 
tion. An interesting possibility that we shall soon investigate 
with multi-mass models is the concentration at the centre 
of more massive stars, i.e. mass segregation. 

Note that some anisotropy has been detected among the 
stars orbiting the central massive black hole of the Milky 
Way, Sgr A* , at distances closer than 1 , i.e. 0.04 pc which i s 
well inside the critical radius (> lpc) llSchodel et al.ll2003l) . 
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Figure 7. Profile of the three-dimensional velocity dispersion for 
the same case and time as Fig.|H] See text for comments. 
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Figure 8. Profile of the anisotropy parameter for the same case 
and time as Fig.|S| The decrease of at the border is an artefact of 
the inappropriate boundary condition. An outer boundary with 
radial anisotropy should be open, but here we enforce the adia- 
batic wall. If one "opens" the wall, but it would be at the expense 
of the stability of the program. All this does affect only a very 
small fraction of the total mass. 
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Figure 9. Projected density for the same case and time as Fig.l?fl 
In the interval between r n and r cr ; t we get a slope of —3/4, as 
expected. 
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Figure 10. Projected velocity dispersions for the same case and 
time as Fig.|H] The line of sight component is represented with a 
solid line and the proper motion component with a dashed one 
for the radial and dashed-dot for the tangential contribution. 



However, the detected anisotropy is in the radial direction 
rather than tangential. It is probably not connected to loss- 
cone effects but to particular h istory of these seemingly very 
young stars (iGhez et alJl20odT) which remains a puzzle. 

In Fig. ll3l the diffusion model is layed out for the loss- 
cone as evaluated in previous sections. We evaluate a model 
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Figure 11. Evolution of the stellar density in the central region 
for our model with 10 5 stars and A^bh(0) = 50 Mq (case a). 
The solid line depicts the density at the influence radius R n . The 
dashed line shows the average density within i? n . 
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Figure 12. Same as Fig. llll but for the three dimensional velocity 
dispersion. 



close to the post-collapse moment analysed in the other 
plots. We depict here the loss-cone filling factor K (upper 
panel), the loss-cone and diffusion angles 6o and 9\ c (middle 
panel) and the local contributions to the total loss-cone star 
accretion rate (lo wer panel). Our diffu sion model reproduce 
well the picture of lFrank fc Reesl lll976ft : the critical radius is 



defined by 9d = #i c and coincides with the radius where the 
local contribution to the loss-cone accretion rate has its peak 
value. These two angles are connected to the time-scales t ln 
and t ou t, 0% oc t out /t Icl!ix and 8? c oc imArciax- The figures 
show that the maximum contribution to the mass accre- 
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tion rate stems at th e radius where 8d = #ic- Consistently, 
iFrank fc Reesl lll976l) estimated the total mass accretion rate 

as .Mbh oc p(r c rit ) ^crit / ^rclax (rcrit ) • 

The dependence of loss-cone accretion rate on time 
during the late re-expansion phase can be estimated 
th rough simple scaling laws. One starts with the relation 
of lFrank fc Reesl l)l97ftl mentioned above, 



Mi 



£ relax (?"crit ) 

and the definition of the critical radius, 

^(Vcrit) = #D(Vcrit) 

ft ^ ^cross (?~crit ) 
fait ^relax (^*crit) 

One substitutes the following relations into Eq. 1461 



(45) 



(46) 



trclax(r) OC 



„3/2 



M 



1/2 ' 



(47) 



M 



3/2 



rt(r) r 3 / 2 n(r) ' 



where we have made use of the fact that the potential 
is dominated by the BH in the region of interest (r cr it < rh). 
Finally one needs the dependence of the density of stars on 
time and radius, n(r,T). We have seen that, to a good ap- 
proximation, the re-expansion of the cluster is homologous 
with Lagrange radii expanding like R oc T 2//3 . In the re- 
gion between r cr i t and rh, the density profile resembles a 
power-law cusp; hence, a general self-similar evolution can 
be described by 



T = 9.95x10' 7 rh (0) = 1. 04x10'° yrs 
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Figure 13. In this triple diagram we evince the dependence on 
radius of the mass accretion rate, the loss-cone and diffusion angle 
which are related to the filling and depletion time-scales of the 
loss-cone (see text) and the filling factor K. The critical radius is 
defined by the condition Bp = d\ c (broken line). 



n{r,T) = n {T) 



(48) 



where ro is some Lagrange radius. Hence, from conservation 
of mass inside ro, 



n{r,T) oc T^~r~ a . 
Combining relations 1461 1471 and 1491 one finds 



rent oc M££- 
and, inserting this into Eg. 1451 



.Mbh oc M 
which, by integration, yields 



6(4-c) 
bh 
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» 7 1 3(4-or) 



7(a-3) 
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M bh ocTT3^- ; 



M 



bhtxT^ 



(49) 



(50) 



(51) 



(52) 
(53) 



For a = 7/4, which is appropriate here (again because 
r C rit < rh), the exponent in the last relation turns out to be 
—95/79 ~ —1.20, in remarkable agreement with figures |5Jd 
and|3J>. 



5 DISCUSSION 

We have presented in this work a method to follow the 
evolution of a spherical stellar cluster with a central ac- 
creting BH in a fully self-consistent manner concerning 
the spatial resolution. As regards the velo city space, we 
use a simplified model based on ideas of IFrank fc Reesl 
lll976h in order to describe the behaviour of the distribu- 
tion function inside and outside the loss-cone by a simple 
diffusion equation. This numerical method is an extension 
of the "gaseous model" which has been successfully ap- 
plied to a variety of aspects of the evolution of globular 
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ICiersz fc Spurzeml 1200(1 12003 



version, the simulation of galactic nuclei is also feasible. 

In addition to an explanation of the physical and nu- 
merical principles underlying our approach, we have concen- 
trated on a few simple test computations, aimed at check- 
ing the proper behaviour of the code. We considered a sys- 
tem where all stars are and remain single, have the same 
mass, stellar evolution and collisions are neglected and a 
seed central BH is allowed to grow by accreting stellar mat- 
ter through tidal disruptions. The present version of the 
code already allows for a (discretised) stellar mass spectrum 
and stellar evolution and we are in the process of includ- 
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ing stellar collisions because they are thought to dominate 
over tidal disruption in most galacti c nuclei, as far as a c- 
cretion on to the BH is concerned jDayid et all Il987al lbt 
iMurphv et~al1ll99ll iFreitag fc Bendl2002l) . In a subsequent 
paper, we shall increase complexity and realism one step 
further and consider systems with a mass spectrum. Us- 
ing both this gaseous c ode a nd the Monte Carlo algorithm 
jFreitag fc Bend 1200 ll . 120021) . we will investigate the role 
of mass segregation around a massive black hole (Amaro- 
Seoane, Freitag & Spurzem, in preparation), a mechanism 
which may have important observational consequences as 
it probably affects the structu r e of the central cluster 
of th e Milky Way dMorrisl Il993l: fMiralda-Escude fc Gouldl 
120001: iFreitad l2003blat IPfahl fc Loebl 12003ft and impacts 
rates of tidal disruptions and capture of compact stars 
by em ission of gravitational wave s in dense galactic 
nuclei l|Ma gorrian fc Tremainel Il999l: ISver fc lllmed" Il999t 
ISigurdssonll2003l . and references therein). 

Unfortunately, the literature has relatively little to offer 
to check our models. The most robust predictions are proba- 
bly the analytical and semi-analytical analysis for the regime 
where the gravity of the BH dominates and the Fokker- 
Planck treatment of relaxation holds (lBalK^Jl_fc_Wolj 
197rllShapfro fc Lightmanlll97rllLightman fc Shaoir Jl977t 
Cohn fc Kulsnidll978ri . The most important feature of these 
solutions is that, provided the system is well relaxed and one 
stands beyond the critical radius (inside of which loss-cone 
effects complicate the picture), a cuspy density distribution 
is established, p oc r~ a with a = 7/4. Our code nicely agrees 
with this prediction. 

Concerning the evolution of the system, we first note 
that, initially, the cluster follows the usual and well un- 
derstood route to core-collapse. That the gaseous model 
can successfully simulate this phase has been clearly 
established in previous works iGiersz fc Spurzeml Il994t 
ISpurzem fc Aarsethll996l) . When the core has become dense 
enough, the BH starts growing quite suddenly. As it ac- 
cretes stars that are deeply bound, i.e. with very negative 
energies, the BH creates a outward flux of energy and al- 
lows the cluster to re-expand. As long as the source of 
energy is centrally concentrated and that the mass of the 
BH remains relatively slow, one expect the re-expansion 
to become self-similar, a regime dur ing which the size of 
the cluster increases like R oc T 2/3 faenonlll965t Ishapirol 
ll977l:lMcMillan et al.iri98ll : lGoodmanlll984l . among others). 
This is again well reproduced by the gaseous model. Solv- 
ing the Fokker-Planck equat ion w ith a Monte Carlo method , 
IMarchant fc Shapird lll980t) and iDuncan fc Shapird lll982T) 
have realised a series of simulations of single-mass globu- 
lar clusters with a central BH. Because their resolution was 
quite low and because they used "initial" conditions diffi- 
cult to implement (in most of their runs the central BH is 
not present initially but introduced at some instant dur- 
ing deep collapse), we do not attempt a quantitative com- 
parison with their results. An added difficulty is that we 
do not include tidal truncation of the clus t er. Ho wever, an 
important finding of IMarchant fc Shap iro ( 198Cj) is repro- 
duced by our computations, namely that the initial mass 
of the seed black hole has little effect on the post-collapse 
evolution, provided it represents only a small fraction of 
cluster mass. In particular, the BH mass at late times con- 
verges to the same value which only depends on the size 



and mass of the cluster. We note that such convergence was 
also obtained with the Monte Carlo algorithm and that a 
comparison between results obtained with that code and an 
early version of t he pro gram described here was presented in 
IFreitag fc Bend J2002f) . More comparisons between the two 
methods are planned (Amaro-Seoane, Freitag & Spurzem, 
in preparation). 

Here one has to mention that the energy input of the 
BH star accretion causes a "temperature" increase in the 
central region which is followed by a thermal expansion. 
Therefore, the system is a normal thermal system with pos- 
itive specific heat in contrast to the cores of self-gravitating 
systems, where the energy input due to binary harden- 
ing causes a core expansion and decrease of temperature 
iBettwieser fc Sugimotd lT984r). Afterwards, an inversion in 
the radial temperature profile foll ows and the ex pansion is a 
reverse gravothermal instability iSpurzemlll99 if) . Since our 
system is dominated by external gravitation in the centre we 
cannot expect such a behaviour. Furthermore, as the central 
BH grows irreversibly, it continues accreting stars in spite 
of the re-expansion and continuous decrease of the density 
of stars. Hence, a second core collapse is impossible and os- 
cillations of the central cluster density do not occur. 

Among the aspects of our results that require further 
investigation, we mention the shape of the density profile in- 
side the critical radius. Although the well known p oc r~ 7//4 
Bahcall-Wolf solution only strictly applies for r cr i t < r < rh, 
this cusp extends inward nearly dow n to the tidal disrup- 
tion r adius in the stationary models of lMarchant fc Shapird 
il979F) . which also include loss-cone physics. This result is 
in dis agreemen t with the analysis of [Pokuchaev fc Ozernoil 
Jl977ft (see also[5 zernoi fc Reinhardtll978r) which predicted 
p oc r _1//2 in this region. As shown on Fig. |S| we obtain an 
even stronger flattening of the density law inside r C rit- At the 
present time it is unknown to us which solution, if any, is the 
correct one. A possibility to be considered is that this a con- 
sequence of the truncation of the moment equations to the 
second order. In other words, in regions were the loss-cone 
is significantly depleted, representing the velocity distribu- 
tion by a simple dispersion ellipsoid and using the veloc- 
ity dispersion to determine an "effective" loss-cone aperture 
(Eq. [SJ is clearly quite a strong approximation. This may 
impact the density distribution as the system adjusts its 
central structure to produce the heating rate required by 
the overall expansion. 

Fortunately, numerical inaccuracies at very small radii 
are unlikely to affect the overall structure and evolution of 
the cluster because the loss-cone accretion physics are es- 
sentially determined by the conditions at the critical radius, 
not in the immediate vicinity of the BH. 



APPENDIX A: A BRIEF MATHEMATICAL 
DESCRIPTION OF THE CODE 

In this appendix, we explain briefly how the gaseous model 
is solved numerically. We concentrate on aspects of the 
method not exposed in previous paper s. This description 
is ther efore complementary to Sec. 2.2 of lGiersz fc Spurzeml 
ill 9941) . The algorithm used is a pa rtially implicit New ton- 
Rap hson-Henyey iterative scheme (iHenvev et al.l fi~959l . see 
also lKippenhahn fc Weigert1ll994l. Sec. II. 2). 
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Figure Al. Representation of the logarithmic radial mesh used 
in the code. With v, p we represent both, the radial and tangential 
components of the velocity and pressure. 



Putting aside the bounding conditions, the set of equa- 
tions to be solved are Eqs. 1351 to 1381 In practice, how- 
ever, the equations are rewritten using the logarithm of all 
positive quantities as dep endant functions. As explained in 
lOiersz fc Spurzernl l|l994l ) , this greatly improves energy con- 
servation. Formally, one may write this system as follows 



dt 



+ f 
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f (0 
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dx (j) ] 
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= for i = 1 . 



= for i = 5. 
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(Al) 

where the x^ are the local quantities defining the state of 
the cluster, i.e. 



x = {x m ,x^,...x (N ^} 

= {log P, U, log p r , log p t , log M r 



(A2) 



To be solved numerically, this set of coupled partial dif- 
ferential equations have to be discretised according to time 
and radius. Let us first consider time stepping. Let At be the 
time step. Assume we know the solution x(t — At) at time 
t — At and want to compute x(t). For the sake of numerical 
stability, a partially implicit scheme is used. We adopt the 
shorthand notations x (i) = x {i) (t) and y {i) = x {i) {t - At). 
Time derivation is replaced by finite differences, 

dx (i) 



dt 



At~ 



(A3) 



In the terms / 



(0 



we replace the x'^ by x^ which are values 



intermediate between and x^\ — Cx^ + (1 
with £ = 0.55 for stability purpose (iGiersz fc Spurzernl 
11994) . 

Spatial discretisation is done by defining all quantities 
(at a given time) on a radial mesh, {ri, r2, . . . fjv r } with r\ = 
and rjv r = r max . A staggered mesh is implemented. While 



values of r, u, vt, v r and M r are defined at the boundaries of 
the mesh cells, p, pt and p r are defined at the centre of each 
cell, see Fig. lAll When the value of a "boundary" quantity 
is needed at the centre of a cell, or vice-versa, one resort 
to simple averaging, i.e. b k = 0.5(6fe_i + b k ), c k = 0.5(c fe + 
Cfc + i), if b and c are is border- and centre-defined quantities, 
and b, c their centre- and border-interpolations, respectively. 
For all runs presented here, N T = 300 and the points ri to 
f'max are logarithmically equidistant with r max = 10 4 pc and 
ri ~ 1.7 x 10 -6 pc. Let us adopt the notation x^' 
value of x at position r k (or f k ) and Ar^ 
Then, radial derivatives in the terms are approximated 
by finite differences, 



»y for the 
r k - r-fe-i. 
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if the derivative has to be evaluated at a point where Xk is 
defined (centre or border of a cell), or 



dx u) 
dr 
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x k 



? 0) 
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otherwise. As an exception we use upstream differencing in 
du/dr for the second equation in set 1351 i.e. the difference 
quotient is displaced by half a mesh point upstream to im- 
prove stability. 

By making the substitutions for dx^ /dt and dx^ /dr 
in the set of differential equations IA1I one obtains, at each 
mesh point r k , a set of jV eq non- linear algebraic equations 
linking the new values to be determined, x 



and x k , to 



the "old" ones, 



and y k , 



which are known, 



•> k 



1 



= 

Nr. 



(A6) 



Note that the structure of the equations is the same at 
all mesh points, except k = 1 and k = N T . In particular, 
terms k — 1 do not appear in . Also, one has to keep in 
mind that only the x_ k -i ano - £fc are unknown; the y k and 
y play the role of fixed parameters in these equation (as do 
the Ar k ). If one defines a (-/V oq x N T )- dimension vector X* 
whose component N cq (k — 1) + i is x k , one can write the 



system of iV eq x iV r equations as T*{X*) = 0, i.e 
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The system is solved iteratively using Newton-Raphson 
scheme. If X*i m ] is the approximation to the solution of 
Eg. IA7l after iteration m, with T*\ m ] = T* (X* r m i) 7^ 0, the 
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solution is refined through the relation 

(dT* 



X [m + l] — X [ 



\8X* 



(A8) 



where (dT* /dX*)^ 1 is the inverse of the matrix of deriva- 
tives. The latter, of dimension (N eq N T ) x (AT eq N T ), has the 
following structure 
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In this diagram, each square is a AT oq x N aq sub-matrix. For 
2 ^ k ^ N T — 1, lines N cq k — 6 to N e<1 k of dT* /dX* are com- 
posed of a group of 3 such iV eq x 7V eq matrices, 0_ k , ■ >,. C+ k 
that span columns N eq k — 13 to N cq k + iV eq , while the rest 
is composed of zeros, 
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The Heyney method is a way to take advantage of the 
special structure of matrix dT*/dX* to solve system TA8I 
efficiently, with a number of operation scaling like 0(N T ) 
rather than O(N^) as would be the case if one uses a general- 
purpose matrix inversion scheme 6 . Setting B* = —T*[ m ] and 
W* = X*[ m+ i] — X*[ m ], Eg. IA8I is equivalent to 



(If 1 



(All) 



with W* the unknown vector. We further decompose vec- 
tors W* and B* into N eq -dimensional sub-vectors, each one 
representing the values at a given mesh points, 



W* 



W 2 



(A12) 



\WnJ 

6 Memory usage is also reduced, scaling like 0(N r ) rather than 
<D(N?). 



Then, the system IA12I can be written as a set of coupled 
Aeq-dimensional vector equations, 

■ 1W1 + a +1 W 2 = Si 

□- fc Wfe-i + m k w k + a +k w k+1 = B k (Ai3) 

jv r Wjv r -i +»n, Wjv r = B Nr ■ 

The algorithm operates in two passes. First, going from k = 
1 to JVi, one defines recursively a sequence of A^q-vectors 
Vfe and (N Bq x Ae^-matrices SfJtfc through 

Vi = (Bi^Bi 

smi = («i)- 1 n +1 

V fe = (U k - □_ fc OTfc_i)" 1 (B k - D-feVfe-i) 



(A14) 
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2 < k < N r 



9Jijv r is not defined. In the second pass, the values of the 
unknown V k are computed, climbing back from k = N T to 
1, with 



Wjv r = V Nr 
W fc = V fc - OT fc W fc+ i 



1 < k < AT r - 1. 



(A15) 



Note that, with this algorithm, only (Neq x A" eq ) matrices 
have to be inverted. We use Gauss elimination for this pur- 
pose because this venerable technique proves to be robust 
enough to properly deal with the kind of badly conditioned 
matrices that often appear in this application. 

The initial model for the Newton- Raphson algorithm is 
given by the structure of the cluster at the previous time, 
X*[o](t) = X*(t—At) One iterates until the following conver- 
gence criteria are met. Let us set Sx 
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-i] 



Then, the condition for logarithmic quantities is 
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with ei = 10 6 . For velocities (u, v r — u, Vt — u), one checks 
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< £2, 
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with £2 = 10 3 and w k = r k (4nGp k ) 1 ^ 2 . Generally, two 
iterations are sufficient to reach convergence. 



APPENDIX B: VELOCITY DISPERSION IN 
THE CENTRAL REGIONS 

In the region dominated by the central BH, one may expect 
a "Keplerian" profile for the velocity dispersion, a oc r _1//2 . 
However, in Sec.|lJ we have seen that in our standard model 
(10 5 stars, case a this relation does not really apply where 
expected, i.e. between the "1-star" and the influence radius. 
Here we show that the spherical Jeans equation for a system 
in dynamical equilibrium is nevertheless obeyed. 

In spherical symmetry the assumption of dynamical 
equilibrium (or stationarity) amounts to u = (v T ) = 0. We 
note that the gaseous model can cope with a / 0. On the 
other hand, for a system whose evolution is driven by relax- 
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ation, one expects a r .t S> u . The Jeans equation then reads 
jBinnev fc Tremaindfl987t Eq. 4-55) 



GM T = -rat 



d In n d In a 2 
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+ 2/3 



(Bl) 



M r is the mass enclosed by the radius r, n the number den- 
sity of stars and 2(3 — 2 — a 2 /a 2 is the anisotropy parameter; 
other quantities have been defined previously. One sees eas- 
ily that if M r = Mhh and the first and third term in the 



brackets are both constant, 



-1/2 



at small r. But, as 



figures |B"T1 and lB"2l demonstrate, none of these assumptions 
exactly applies in the range of radius under consideration. 
Consequently, we do not get a clean Keplerian velocity pro- 
file although (or because) Ea . IB H is satisfied. Finally, we men- 
tion that our models with 10 6 and 10 7 stars (and same initial 
size) exhibit a Keplerian velocity cusp outside the 1-star ra- 
dius, during their post-collapse evolution. This is partly due 
to the relatively more massive black hole (larger influence 
radius) and partly to the much smaller 1-star radius. 



APPENDIX C: MBH WANDERING IN A 
CUSPY CLUSTER. 

Here we present a simple estimate of the wandering radius 
fj wan of a MBH embeded in a stellar cluster whose density 
posseses a power-law cusp in the inner regions. We assume 
that, were it not for the effect of the MBH it self, the stellar 
cluster would be desc ribed by an eta-model JPehnenlll99l 
iTremaine et alJll994l) with enclosed stellar mass 

M, { r)=M cl { J ^)\ (CI) 

where M c \ is the total mass in stars and a the break radius. 
For r <C a, the density is p oc r~ a with a = 3 — r\. Inside 
-Rwan, the MBH strongly perturbates stellar orbits and we 
suppose the density is rendered more or less constant. Hence, 
the potential felt by the MBH is approximately harmonic, 



$(r) = $o + ^ 2 r 2 , with to 2 = 



R 3 

J 1 vv.'l I I 



(C2) 



For a harmonic oscillator, the RMS amplitude of the oscil- 
lations in velocity and space are linked to each other, 



Vrms — w -Rrms 



k-'Rwan — 



GM*(R* 

i?wan 



(C3) 



iDorband et alj (120031) have verified with iV-body simulations 
that equipartition of kinetic energy between the MBH and 
the stars is established, at least in the case of r\ — 1.5. 
Namely, 



(C4) 



whe re a is the stellar velo city dispersion at r = a. For 77 
1.5 jTremaine et alJll994h . 



2 ni GMd 
a ~ 0.1 . 



(C5) 



Finally, assuming i? W an <S a and, hence, X,(E W an) — 
Mcl {Rwan/a) v and combining equations IC3I IC~4l and lUsl 

we obtain 

2 



-Rwan « 0.01a 



for tj = 1.5 (C6) 



and 



Rw 



1/(1-1) 



(C7) 



for general eta-models. 



APPENDIX D: PROJECTED VELOCITY 
DISPERSIONS 

For the Fiaures HCJl andl9lwe have integrated the density along 
the z-axis for the projection, 

/•z=il max 

E(r)= / ,o(V V 2 + z 2 ) dz = 



2=0 

z=R„ 



P (R) 



R 



v^R 2 



■dR 



(Dl) 



If one observes the cluster along the z-axis, the contribu- 
tions to the projected velocity dispersions are as indicated 
in Fig. IDH We have that R = \/r 2 + z 2 and ag — = at, 
where the subscript t stands for tangential. We can reckon 
that 



2 2 2 /1 1 ~ 2 _• 2 /1 

<7 Z = <7r cos a + a t sm a 

°t = a R sm2 6 + 5\ cos 2 8, (D2) 

where we have defined a 2 = a 2 /2, to be consistent with the 
notation used until now. 

In Fig.lUTl we can see that <r z contributes to ctlos, to 
(T pnliI and <t 2 to <T pm ,t. 

Thus, we obtain the projected velocity dispersions, 



S(R) J z= 



f R ™\al{R) 

J z=0 



cos 6 + & 2 (R) sin 6)p(R) dz 



{jp{ai(R)-at{R)) + at{R))p{R)dz 



0"pm,r(r) = 



" (<x R (-R) sin 2 + &t(R) cos 2 0)p{R) dz 



S(R) 



2 /-flmax 

V>) = f^y J a 2 {R)p(R)dz, 

since sin 2 6 = r 2 /R 2 = 1 - cos 2 9 and cos 2 9 = z 2 /R 2 



(D3) 
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Figure Bl. Check of stationary Jeans equation for our "stan- 
dard" model (10 5 stars, case a) at T = 10 Gyrs (same model and 
time as Fig.|7J. Panel (a) depicts the logarithmic derivatives of 
the stellar density n and radial component of the velocity dis- 
persion well as the anisotropy parameter, 2/3 = 2 — cr^/cr^. Hor- 
izontal lines corresponding to values of 0, —1, —1.75 and —2.23 
are present to guide the eye. In particular, the "Keplerian" ve- 
locity profile is din o^/dln r = —1. The round dot indicates the 
influence radius, the star the critical radius and the square the 
"1-star" radius. On panel (b), we plot the three terms of the right 
side of the stationary Jeans equation and check that their sum is 
(nearly) equal to the left side term, i.e. GM T . 
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Figure Dl. Different contributions to the projected velocity dis- 
persions. The semi-major axis of the ellipsoid perpendicular to 
the page corresponds to <rt = crj,. 



